Spatial Analysis
Julia Earle, Fatemeh Ghaheri, and Angelina Locker
May 4, 2019
Preliminaries
Load and install the following packages: {tmap}, {leaflet}, {rgdal}, {rgeos}, {sp}, {raster}, {adehabitatHR}, {devtools}, {rayshader}, {mapmate}
Objectives
The objective of this module is to demonstrate some of the capabilities R has for handling spatial data. We will provide you with a brief overview of the various formats of spatial data, explain how to import spatial data into R, and provide a few examples of spatial data manipulation.
Before we begin with the data, you’ll need to link up the the GitHub repository. Create a New Project –> Version Control –> Git: https://github.com/ajlocker/ADA-SpatialAnalysis.git
This will allow you to access all the necessary files for this module.
Before we begin with the data, you’ll need to link up the the GitHub repository. Create a New Project –> Version Control –> Git: https://github.com/ajlocker/ADA-SpatialAnalysis.git
This will allow you to access all the necessary files for this module.
Shapefile
Shapefile is a data storage format where you can store different types of geographic information and geospatial data such as location, shape, line and points and their various details such as names etc. The information in a shapefile is actually a big set of files with large data in them combined within one shapefile. Like a csv file shapefile (shp file for short) also has different set of details and files where large data is stored that can be seen in columns when loaded to R. Many GIS applications and softwares such as ArcGIS use shp file for their visualizations, mapping and analyses. But since they are not usually free it’s a great feature that R has as a free tool.
Point data vs. raster data
Point data represent distinct and separate points and are commonly defined by their geographic coordination and can be located within a polygon, for example location of a bridge and an archaeological site can be shown as a point.
Raster data or “gridded” data represents surfaces and are saved as pixels where each pixel represents a surface.
So now that we know the definition of the above terms we know that I mentioned them because they are among the most important features in geospatial visualization.
First we will show how you can load GIS data into R, and how you can manipulate it to create basic maps.
library(tmap)## Warning: package 'tmap' was built under R version 3.5.3library(leaflet)## Warning: package 'leaflet' was built under R version 3.5.3library(rgdal)## Warning: package 'rgdal' was built under R version 3.5.3## Loading required package: sp## Warning: package 'sp' was built under R version 3.5.3## rgdal: version: 1.4-3, (SVN revision 828)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/ajlocker/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/ajlocker/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1library(rgeos)## Warning: package 'rgeos' was built under R version 3.5.3## rgeos version: 0.4-3, (SVN revision 595)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-1
## Polygon checking: TRUElibrary(sp)
library(raster)## Warning: package 'raster' was built under R version 3.5.3library(adehabitatHR)## Warning: package 'adehabitatHR' was built under R version 3.5.3## Loading required package: deldir## deldir 0.1-16## Loading required package: ade4## Warning: package 'ade4' was built under R version 3.5.3## Loading required package: adehabitatMA## Warning: package 'adehabitatMA' was built under R version 3.5.3##
## Attaching package: 'adehabitatMA'## The following object is masked from 'package:raster':
##
## buffer## Loading required package: adehabitatLT## Warning: package 'adehabitatLT' was built under R version 3.5.3## Loading required package: CircStats## Warning: package 'CircStats' was built under R version 3.5.3## Loading required package: MASS## Warning: package 'MASS' was built under R version 3.5.3##
## Attaching package: 'MASS'## The following objects are masked from 'package:raster':
##
## area, select## Loading required package: boot## Warning: package 'boot' was built under R version 3.5.3Using spatial data in R
#Using a one line code you can load shapefile into R and read it: #m <- readOGR(".", "Zipped_SPA_Shape_file")
#The other format of code you can use for this purpose is m1 <- readOGR("shp/Zipped_SPA_Shape_file.shp",layer="Zipped_SPA_Shape_file"), but i found this format a little problematic and confusing, and it also gives erros if you dont use it carefully, so i recomment the first code
#The “.” Is the source and the "Zipped_SPA_Shape_file" is the layer which is the name of the file without the suffixes.
#read the shapefile using the bellow code
# Using a one line code you can load shapefile into R and read it: m <- readOGR(".", "Zipped_SPA_Shape_file")
# The other format of code you can use for this purpose is m1 <- readOGR("shp/Zipped_SPA_Shape_file.shp",layer="Zipped_SPA_Shape_file"), but i found this format a little problematic and confusing, and it also gives erros if you dont use it carefully, so i recomment the first code
# The “.” Is the source and the "Zipped_SPA_Shape_file" is the layer which is the name of the file without the suffixes.
# read the shapefile using the bellow code
mydata1 <-readOGR(".", "Zipped_SPA_Shape_file")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\ajlocker\Development\Projects\DataAnalysisClass\repos\ADA-SpatialAnalysis", layer: "Zipped_SPA_Shape_file"
## with 6 features
## It has 15 fields
## Integer64 fields read as strings: FID OBJECTIDmydata1## Object of class "SpatialPolygonsDataFrame" (package sp):
##
## Number of SpatialPolygons: 6
##
## Variables measured:
## FID OBJECTID MAPKEY SQ_MILES ACRES SPA
## 0 1 1 2B 71.19565 45565.22 SPECIAL PLANNING AREA 6
## 1 2 2 2B 42.99561 27517.19 SPECIAL PLANNING AREA 1
## 2 3 3 2B 58.89297 37691.50 SPECIAL PLANNING AREA 3
## 3 4 4 2B 23.98706 15351.72 SPECIAL PLANNING AREA 2
## 4 5 5 2B 57.59823 36862.87 SPECIAL PLANNING AREA 5
## 5 6 6 2B 47.13667 30167.47 SPECIAL PLANNING AREA 4
## GlobalID created_us created_da last_edite
## 0 {3FFECD3C-A2A4-428E-BE34-BC78AF729445} <NA> <NA> <NA>
## 1 {0930519B-24B9-42F5-A509-C985F03A88C4} <NA> <NA> <NA>
## 2 {E28E8320-B988-4BF1-8212-32710A8B21C1} <NA> <NA> <NA>
## 3 {0FEB811D-0B9D-4F4F-AD9B-E5C4BA4A7C77} <NA> <NA> <NA>
## 4 {3A32692E-8248-405B-B428-F78A97FB900B} <NA> <NA> <NA>
## 5 {FCF1AC76-46DD-4EFA-B8D4-5032AAAE017A} <NA> <NA> <NA>
## last_edi_1 Shape_STAr Shape_STLe Shape__Are Shape__Len
## 0 <NA> 1984828850 203611.7 268353008 74823.16
## 1 <NA> 1199698478 182234.3 161187360 66796.64
## 2 <NA> 1640887228 226639.7 220714059 83128.26
## 3 <NA> 668854756 123557.2 90078255 45316.48
## 4 <NA> 1605752865 185957.2 216438440 68304.16
## 5 <NA> 1314100339 191163.3 177292582 70159.31head(mydata1)## FID OBJECTID MAPKEY SQ_MILES ACRES SPA
## 0 1 1 2B 71.19565 45565.22 SPECIAL PLANNING AREA 6
## 1 2 2 2B 42.99561 27517.19 SPECIAL PLANNING AREA 1
## 2 3 3 2B 58.89297 37691.50 SPECIAL PLANNING AREA 3
## 3 4 4 2B 23.98706 15351.72 SPECIAL PLANNING AREA 2
## 4 5 5 2B 57.59823 36862.87 SPECIAL PLANNING AREA 5
## 5 6 6 2B 47.13667 30167.47 SPECIAL PLANNING AREA 4
## GlobalID created_us created_da last_edite
## 0 {3FFECD3C-A2A4-428E-BE34-BC78AF729445} <NA> <NA> <NA>
## 1 {0930519B-24B9-42F5-A509-C985F03A88C4} <NA> <NA> <NA>
## 2 {E28E8320-B988-4BF1-8212-32710A8B21C1} <NA> <NA> <NA>
## 3 {0FEB811D-0B9D-4F4F-AD9B-E5C4BA4A7C77} <NA> <NA> <NA>
## 4 {3A32692E-8248-405B-B428-F78A97FB900B} <NA> <NA> <NA>
## 5 {FCF1AC76-46DD-4EFA-B8D4-5032AAAE017A} <NA> <NA> <NA>
## last_edi_1 Shape_STAr Shape_STLe Shape__Are Shape__Len
## 0 <NA> 1984828850 203611.7 268353008 74823.16
## 1 <NA> 1199698478 182234.3 161187360 66796.64
## 2 <NA> 1640887228 226639.7 220714059 83128.26
## 3 <NA> 668854756 123557.2 90078255 45316.48
## 4 <NA> 1605752865 185957.2 216438440 68304.16
## 5 <NA> 1314100339 191163.3 177292582 70159.31# Plot the data into a map
plot(mydata1)
# Creates a choropleth map of our any( in this case SQ_MILES) variable
# With a choropleth map you can show statistical data with different specifications such as color and symbols shading on the map
## qtm function: allows you draw a thematic map. when you add fill= to the code then you specify how you want your details to be shown on the map. but without fill= you will just get a map without any specification.
m <- qtm(mydata1, fill = "SQ_MILES")
m## Warning: package 'sf' was built under R version 3.5.3## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
# Removes the black borders
# here with tm_shape fucntion we specify a spatial data object and then we can add various elements to it.
tm_shape(mydata1) + tm_fill("SQ_MILES")
# Add borders
# Alpha is the transparency of the borders
# using the function tm_borders we are able to add borders.
tm_shape(mydata1) + tm_fill("SQ_MILES") + tm_borders(alpha=.9)
# Adding a north arrow
tm_shape(mydata1) + tm_fill("SQ_MILES") + tm_compass()
# Giving a title to the map
tm_shape(mydata1) + tm_fill("SQ_MILES", palette = "Reds",
style = "quantile", title = "% with a SQ_MILES") +
tm_borders(alpha=.4) +
tm_compass() +
tm_layout(title = "Camden, Arizona", legend.text.size = 1.1,
legend.title.size = 1.4, legend.position = c("right", "top"), frame = FALSE)## Legend labels were too wide. The labels have been resized to 0.97, 0.97, 0.97, 0.97, 0.97. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Density Mapping
In the density mapping section, we will learn how to visualize concentrations of point data. Specifically, we will use kernel density estimation to calculate where the highest concentrations of point data are located in our data set.
Density mapping is a means of visualizing concentrations of spatial phenomena. Point data are most often used for density mapping, but lines can also be used. If you have a high volume of point data distributed across a large area, it may be difficult to understand the distribution of points. Density maps can simplify this data by showing where points are concentrated.
Density mapping can be incredibly useful when you have a large data set representing isolated phenomena over a continuous surface (e.g., artifacts, populations). As such, density mapping would not be useful for representing small data sets, or data that represents continuous phenomena (e.g., elevation or temperature). With spatial data visualization it is important to understand first what kind of data you have, and what you want to represent or analyze.
# Download data from https://github.com/ajlocker/ADA-SpatialAnalysis/blob/master/DEPARTAMENTO.shp
# Remember, the files must be placed into your working directory
# Loading and plotting shapefile (Peru and its departamentos)
DEPARTAMENTO <- readOGR(".", "DEPARTAMENTO")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\ajlocker\Development\Projects\DataAnalysisClass\repos\ADA-SpatialAnalysis", layer: "DEPARTAMENTO"
## with 24 features
## It has 3 fields# Setting the coordinate system for the shapefile
proj4string(DEPARTAMENTO) <- CRS("+init=EPSG:32718")## Warning in `proj4string<-`(`*tmp*`, value = new("CRS", projargs = "+init=EPSG:32718 +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransformplot(DEPARTAMENTO)
This shapefile represents a map of Peru with the departments outlined. We will plot the distribution of recorded archaeological sites in Peru onto this map.
# Download data from https://github.com/ajlocker/ADA-SpatialAnalysis/blob/master/SitioArqueologico.shp
# Loading and plotting point data (recorded archaeological sites in Peru)
SitioArqueologico <- readOGR(".", "SitioArqueologico")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\ajlocker\Development\Projects\DataAnalysisClass\repos\ADA-SpatialAnalysis", layer: "SitioArqueologico"
## with 7907 features
## It has 7 fields# Setting the coordinate system for the point data
proj4string(SitioArqueologico) <- CRS("+init=EPSG:32718")## Warning in `proj4string<-`(`*tmp*`, value = new("CRS", projargs = "+init=EPSG:32718 +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransformplot(SitioArqueologico)
Kernel density estimation
To begin, we will perform kernel density estimation on the point data. Kernel density estimation uses a moving quadrant to calculate the density of points for each area within a given threshold.
For kernel density estimation, be sure you have loaded the package “adehabitatHR” into your library. Here we use the kernelUD function. UD refers to Utilization Distribution, which is a bivariate function. Other functions are available to calculate kernel density in this package.
These points represent the distribution of recorded archaeological sites in Peru. Here we see that we have a very high density of points and it is difficult to distinguish where the higher concentrations of points are located at this resolution.
kde.output <- kernelUD(SitioArqueologico, h="href", grid = 1000)## Warning in kernelUD(SitioArqueologico, h = "href", grid = 1000): xy should contain only one column (the id of the animals)
## id ignored# Grid determines how many pixels will be used for the map. Using a smaller number of pixels will run faster, but will not be as detailed. Grid value should be selected according to how much detail you need to represent.
plot(kde.output)
If we want to map this density estimate, we need to first convert it to a raster.
# Convert to raster
kde <- raster(kde.output)
# Map raster
tm_shape(kde) + tm_raster("ud")
The value in the legend refers to the value calculated from the UD function. A higher UD corresponds to a higher density here.
These maps are much more useful than the original point data to tell us about where most archaeological sites have been recorded. However, without the boundaries of the country of Peru and its departments, it is difficult to determine where exactly the highest densities of archaeological sites are located. To this end, we will plot our kernel density estimates on the original shapefile.
We cannot map the kernel density directly onto the shapefiles, so we must create ‘catchment areas’ at 75%, 50%, and 25% intervals.
# Mask the raster by the output area polygon
masked_kde <- mask(kde, SitioArqueologico)
# Catchment boundaries for kernel density estimates
# Top 75% highest density
range75 <- getverticeshr(kde.output, percent = 75)
plot(range75)
# Top 50% highest density
range50 <- getverticeshr(kde.output, percent = 50)
plot(range50)
# Top 25% highest density
range25 <- getverticeshr(kde.output, percent = 25)
plot(range25)
With this code we have generated our catchment areas, which we can map onto the shape file and point data. Next, we will map the point data and the catchment boundaries onto our original map of Peru with its departments.
tm_shape(DEPARTAMENTO) + tm_fill(col = "#f0f0f0") + tm_borders(alpha=.8, col = "white") + tm_shape(SitioArqueologico) + tm_dots(col = "blue") +
tm_shape(range75) + tm_borders(alpha=.7, col = "#fb6a4a", lwd = 2) + tm_fill(alpha=.1, col = "#fb6a4a") +
tm_shape(range50) + tm_borders(alpha=.7, col = "#de2d26", lwd = 2) +
tm_fill(alpha=.1, col = "#de2d26") +
tm_shape(range25) + tm_borders(alpha=.7, col = "#a50f15", lwd = 2) +
tm_fill(alpha=.1, col = "#a50f15") + tm_layout(frame = FALSE)
To review, with the code outlined above we have calculated kernel density estimates and mapped these densities onto shapefiles and point data in order to understand where the highest densities of points are located on our map.
LiDAR Data and R
Light Detection And Ranging (LiDAR) is very useful survey method that uses pulses of laser light to measure reflected distance between a given sensor (typically on the bottom of an airplane) and a survey area.
LiDAR Illustration
This becomes exceptionally useful when your research area is in the middle of a neo-tropical jungle. 3D images and maps can be created based on the differences in response time, where objects that are higher in elevation are returned more quickly than objects lower in elevation.
Caana Temple, Caracol
You can see that the forest canopy is fairly dense around the site. LiDAR essentially allows us to complete an aerial survey and remove the canopy layer to see what is happening below the tree line. We then have the ability to “groundtruth” the LiDAR data to see if the data match what is actually on the surface.
Caracol map compared to LiDAR imagery
Using LiDAR data to create Digital Elevation Models in R
The first thing we will need to do is set up a working directory for the files.
# Load needed packages
library(raster)
# set working directory to data folder
#setwd("C:/Users/ajlocker/Desktop/R Geospatial")You’ll need to go to download the following files and save them inside your working directory.
https://github.com/ajlocker/ADA-SpatialAnalysis/blob/master/ATXLiDAR/demTIN.tif
https://github.com/ajlocker/ADA-SpatialAnalysis/blob/master/ATXLiDAR/hsTIN.tif
https://github.com/ajlocker/ADA-SpatialAnalysis/blob/master/ATXLiDAR/demlocal.tif
Next, we are going to load in the data by assigning them to be rasters. The DEM file is a Digital Elevation Model, which tells R elevation values from the scanned LiDAR data. The hs file is a Hillslope image.
# assign raster to object
dem <- raster("C:/Users/ajlocker/Desktop/R Geospatial/ATX/demTIN.tif")
# view info about the raster.
dem## class : RasterLayer
## dimensions : 850, 759, 645150 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 620669, 621428, 3348680, 3349530 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=14 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## data source : C:/Users/ajlocker/Desktop/R Geospatial/ATX/demTIN.tif
## names : demTINPlot the Digital Elevation Model
plot(dem, main="Digital Elevation Model \n Austin, Texas")
So here, we can see basic elevation differences from the surveyed area, but it doesn’t show us much detail from the LiDAR imaging. We can examine this data much more effectively by plotting it on top of the Hillshade file, which will essentially add depth to the map.
hs <- raster("C:/Users/ajlocker/Desktop/R Geospatial/ATX/hsTIN.tif")
# view info about the raster.Pay attention to the extent portion. We'll be manipulating this data in a bit.
hs## class : RasterLayer
## dimensions : 850, 759, 645150 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 620669, 621428, 3348680, 3349530 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : C:/Users/ajlocker/Desktop/R Geospatial/ATX/hsTIN.tif
## names : hsTIN
## values : 0, 255 (min, max)# plot hillshade file
plot(hs,
col=grey(1:100/100), # creates a gray scale coloring for the hillshade data
legend=FALSE, # removes any legend
main="LiDAR Survey Austin, Tx", #sets the title of the image
axes=TRUE) # keeps axes titles so we can see the coresponding extent data (xmin, xmax, ymin, ymax)
#then plot the DEM ontop of the hillshade
plot(dem,
axes=T,
alpha=0.5, # transparency of the object (0=transparent, 1=not transparent)
add=T) # add=TRUE (or T); this will add the chm plot on top of the plot we coded immediately above
Creating an Interactive Map
We’re going to use the {leaflet} package to plot our lidar data onto an existing map. We can then zoom in and out to see where our data are geographically.
library(leaflet)
hs <- raster("C:/Users/ajlocker/Desktop/R Geospatial/ATX/hsTIN.tif")
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(hs),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(hs, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(hs),
title = "Elevation")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA# To install the latest version from Github:
library(devtools)## Warning: package 'devtools' was built under R version 3.5.3## Warning: package 'usethis' was built under R version 3.5.3devtools::install_github("tylermorganwall/rayshader")## Skipping install of 'rayshader' from a github remote, the SHA1 (10f48f8b) has not changed since last install.
## Use `force = TRUE` to force installationWe’re going to use a raster file of just elevation points to generate this map. First, I’m going to have you crop the data to make the map a bit more dramatic looking. We’re going to crop the data to a specific quandrant. To do this, we’ll use the {raster} package.
If we look back at our DEM maps, we are going to target the bottom left corner of the image. The axes on the DEM correlate to the extent data from the raster information we pulled. We are going to crop our raster, using the as(extent()) function from the {raster} package
cp <- as(extent(620669, 620900, 3348680, 3349000), 'SpatialPolygons') #Here we are telling R what our new map coordinates will be on a UTM coordinate system and that we want to turn our Raster into spatialpolygons.
crs(cp) <- crs(dem) #crs stands for Coordinate Reference System, and gives us map projection information
atxcp <- crop(dem, cp) #Here we crop the map to be our selected area.
atxcp## class : RasterLayer
## dimensions : 320, 231, 73920 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 620669, 620900, 3348680, 3349000 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=14 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## data source : in memory
## names : demTIN
## values : 62.34, 388.77 (min, max)plot(atxcp) #plot to see the new area. 
Now that our raster is cropped, we will create our 3D maps.
library(rayshader)
#First we need to convert our cropped raster file to a matrix:
atxmat = matrix(raster::extract(atxcp,raster::extent(atxcp),buffer=1000),
nrow=ncol(atxcp),ncol=nrow(atxcp))
#Then we are going to plot with a slight shadow
atxmat %>%
sphere_shade(texture = "bw") %>%
add_shadow(ray_shade(atxmat,lambert = TRUE), 0.7) %>%
plot_map()
#adds ambient shading to the map
ambmat = ambient_shade(atxmat)
#This will generate a 3D interactive map. It should generate a pop-up window that you can then use your mouse to move around as you wish. You can view from the side to see the elevations, zoom in or out, rotate around, etc.
atxmat %>%
sphere_shade(texture = "bw") %>%
add_water(detect_water(atxmat), color="bw") %>%
add_shadow(ray_shade(atxmat,zscale=3,maxsearch = 300),0.5) %>%
add_shadow(ambmat,0.5) %>%
plot_3d(atxmat,zscale=10,fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))
render_snapshot()
You could then take this 3D generated plot and print it on a 3D printer using the command save_3dprint(“atx.stl”)
Concept Review
Spatial Data has sevaral formats and each of them can be read in R.
R has the ability to read and manipulate LiDAR data.
You have the ability to generate interative and 3D maps/models with packages like {leaflet} and {rayshader}.
References
Rayshader.com
https://www.neonscience.org/neon-lidar-rasters-workshop-20150514
https://gis.stackexchange.com/questions/229356/crop-a-raster-file-in-r
Lansley, G., & Cheshire, J. (2016). “An Introduction to Spatial Data Analysis and Visualisation in R”. Consumer Data Research Centre.